Effective Random Matrix Theory description of chaotic Andreev billiards 
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An effective random matrix theory description is developed for the universal gap fluctuations 
and the ensemble averaged density of states of chaotic Andreev billiards for finite Ehrenfest time. 
It yields a very good agreement with the numerical calculation for Sinai-Andreev billiards. A 
systematic linear decrease of the mean field gap with increasing Ehrenfest time te is observed but 
its derivative with respect to te is in between two competing theoretical predictions and close to 
that of the recent numerical calculations for Andreev map. The exponential tail of the density of 
states is interpreted semi-classically. 
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Recently, mesoscopic ballistic two dimensional normal 
(N) dots in contact with a superconductor(S) have been 
extensively studied. Such hybrid systems are commonly 
called Andreev biUiards(li,t2n^i^i5||- the most recent 
works, interest has shifted to mesoscopic fluctuations of 
the excitation spectrum of these systems 9]. Since 

the sub-gap spectrum determines the tunneling conduc- 
tance of an N-S contact this is an essential question both 
experimentally and theoretically. 

Based on the semi-classical treatments and random 
matrix theory (RMT), it was shown by Melsen et al. .2j 
that integrable Andreev billiards are gapless, whereas 
systems with classically chaotic dots possess an energy 
gap on the scale of the Thouless energy Et = ?i/(2r£i), 
where tjj = ttA/{Wv-p) is the mean dwell time in the 
normal dot (here A is the area of the normal dot, W is 
the width of the superconducting region and v-p is the 
Fermi velocity). For such systems, it is assumed that 
5n « Et << a, where = 2nh'^ / (mA) is the mean 
level spacing of the isolated normal dot with effective 
mass TO of the electrons and A is the bulk order param- 
eter of the superconductor 4] . In further studies 6, 7] it 
was concluded that in chaotic cases, the lowest energy 
level El of the system varies from sample to sample with 
a universal probability distribution P{x) given in Ref. 
if the energy levels Ei are rescaled as a; = {Ei — Eg)/ Ag, 
where the mean-field value of the gap Eg and the width 
of the distribution Ag are given by 



Eg = S^"^ = 2-f'l'^ET, 
A„ = Ar" - c'Af i/^Jat 



(la) 
(lb) 



Here 



7 = 5(\/5 — 1) is the golden ratio , c' = 

[(15-6^5) /20]^^V27r, M = Int [kpW/Tr] is the num- 
ber of open channels in the S region and fcp is the Fermi 
wave number (Int[.] stands for the integer part). 

Equations 1^ are strictly valid only in the RMT limit, 
i.e., when the Ehrenfest time te — (1/A) ln(L/AF) tends 
to zero (te is the time needed for a wave packet of mini- 
mal size Af — 27: /kY to spread to the characteristic size L 



of the classically chaotic normal dot with Lyapunov expo- 
nent A). For finite but small enough Ehrenfest time, Sil- 
vestrov et al. 0, and Vavilov and Larking predicted 
that to lowest order in t^/td the mean- field gap Eg de- 
creases linearly by increasing the ratio te/t/). The first 
numerical evidence for the distribution P{x) in the RMT 
limit and the dependence of Eg on the ratio te/tu was 
presented by Jacquod et al. modeling the hybrid sys- 
tem with the one dimensional Andreev map. 

From an experimental point of view, more realistic can- 
didates for studying quantum chaos in hybrid systems 
would be two dimensional Andreev billiards with classi- 
cally chaotic normal region. However, to date no numcri- 




FIG. 1; A normal dot (N) of Sinai billiard in contact with a 
superconductor (S). 

cal study of these systems confirms the predictions of the 
RMT, except for the mean field density of statesj^ . The 
purpose of this paper is twofold. On the one hand, we 
present (for the first time) a numerical study of the gap 
fluctuation and the density of states (DOS) in two dimen- 
sional chaotic Andreev billiards. On the other hand, to 
show that for flnite Ehrenfest time the gap distribution 
P{x) is still the same universal function given in the RMT 
limit but the scaling parameters Eg and A^ extracted 
from the ensemble (see below) are different from those 
given in Eq. (Q. In this way, this is an effective RMT 
description of the fluctuation of the gap using 'renormal- 
ized' scaling parameters. As a further support of our 
model it is also shown that this effective RMT descrip- 
tion correctly predicts the edge of the ensemble averaged 
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DOS. For clarity, we would like to mention that the re- 
cently developed theory [13 for the mean field gap based 
on the concept of a reduced phase space can also be con- 
sidered as an effective RMT description. However, our 
effective RMT description deals with the fluctuation of 
the gap in Andreev billiards. 

In our numerical work we calculated the exact en- 
ergy levels of the so-called Sinai-Andreev (SA) billiard 
in which the normal dot is a Sinai billiard (see Fig. 
The energy levels of the Andreev billiards are the positive 
eigenvalues E (measured from the Fermi energy) of the 
Bogoliubov-de Gennes equation^. To obtain the exact 
energy levels of an Andreev billiard we used the recently 
derived general and quantum mechanically exact secu- 
lar equation expressed in terms of the scattering matrix 
Sq{E) of the normal region^^- The scattering matrix 
So{E) was calculated by expanding the wave function in 
the N region in terms of Bessel functions. 

To ensure that the long classical trajectories (compared 
to the characteristic length of the system) starting and 
ending at the N-S interface are truly chaotic the following 
geometrical constrains should be applied: h = a — W (the 
superconductor is placed at the top of the vertical border 
of the Sinai billiard), R + W >a and R > a/V2. Other- 
wise, there may exist arbitrary long trajectories without 
bouncing on the circular part of the Sinai billiard. For 
such intermittent trajectories the return probability de- 
cays as Pr{s) ~ (s is the length of the trajectory) 
and this results in a gapless energy spectrum of the An- 
dreev bilhardllli. 

For realistic systems the Ehrenfest time is always fi- 
nite and an important question arises, namely how the 
predictions of the RMT for the gap distribution and the 
scaling parameters Eg and Ag are affected. Provided 
that for small enough te/td the distribution of the gap 
in the rescaled variable x remains the same as that in the 
RMT limit, the following equations hold 



Eg = {E,)-{, 

Ag = 5Ei/5x, 



)A 



(2a) 
(2b) 



where {Ei) is the mean value, 5Ei ~ 
is the standard deviation of Ei, while 
5x = \fXx^ 



ViE't) - {Eiy 

x) « 1.21 and 
1.27 are calculated from the gap 
distribution P(a;)[6j. Although the distribution of Ei is 
unknown, its mean and standard deviation can be nu- 
merically estimated from the data of an ensemble of the 
SA billiard, hence Eg and Ag follow. 

Figure 13 shows our numerical results for the integrated 
distribution F[x) = P{x') dx' together with the theo- 
retical prediction (the distribution P{x) is shown in the 
inset). In our numerics we used 5000 slightly different 
realizations of the SA billiard by varying the geometri- 
cal parameters i?, W and the Fermi wave numbers kp 
(for the parameters of the SA billiard see ^3). From 
Eq. Q we found that Eg « 0.5 £'t and Ag « 0.118 £'t 



and they are different from those given in the RMT 
limit, i;™^ = 0.6 £;r and A^"^ = Omi Et- It is clear 
from the figure that the numerical result for F{x) (us- 
ing Eq. O) is also different from that of the theoretical 
prediction in the RMT limit. However, the agreement 
between our numerical results obtained using Eq. Q 
and the universal distribution function F{x) is excellent, 
without any adjustable parameters. Similar numerical 
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FIG. 2: Comparison of our numerical result using the param- 
eters Eg and Ag given by Eq. Q (short dashes) and Eq. (0 
(long dashes) with the theoretical predictions of the integrated 
distribution F{x) (solid line). In the inset the same as in the 
main frame is plotted for the distribution P{x) with square 
and triangular symbols, and solid line, respectively. 

results were found for other ensembles of the SA bil- 
liard. These results imply that for systems with non-zero 
Ehrenfest time the gap distribution is still given by the 
universal function P{x) provided the 'renormalized' pa- 
rameters Eg and A^ obtained from Eq. (j^J are used. One 
can regard the distribution P{x) (and the integrated dis- 
tribution F{x)) with the renormalized parameters as a 
result of an effective RMT description. 

The necessity of the renormalization of Eg and Ag can 
be interpreted as their dependence on the Ehrenfest time. 
Since the characteristic length of the system is uncertain, 
two definitions of the Ehrenfest time, proposed in Ref.llll 
and used in numerical simulations of Ref. were here 
adapted for numerical calculations: 



.(1) 



.(2) 



1 

2A XpLc 
2A^" a7' 



(3a) 
(3b) 



where Lc is the average length of the part of the trajec- 
tory lying between two consecutive bounces at the curved 
boundary segment of the Sinai billiard and Lav = ttA/K 
with perimeter K of the billiard is the mean chord length 
in the normal region. They are parametrically different, 
but their numerical values are of the same magnitude. 
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For the numerical results shown in Fig.|2]the ratio of the 
Ehrenfest time and the dwell time is about t^'' /tjj w 0.1 

(2) 

and ' /m ^ 0-26 (for calculation of the Lyapunov ex- 
ponent see 14j). 

It was predicted theoretically'l3,'ll'| and demonstrated 
numerically using the Andreev map 8] that to lowest or- 
der in t^/td the mean field gap Eg decreases linearly by 
increasing the ratio of the Ehrenfest time and the dwell 
time for ratio much less than one. From our numerical 



0.95 - □ 




0.8 - 

0.75 I • ' • ' • ' • ' • ' • 

0.04 0.06 0.08 0.1 0.12 0.14 0.16 

FIG. 3: The mean field gap Eg/E^'^'^ as a function of ratio 
/rn- The straight line is the numerical fit using Eq. 

study a systematic decrease of Eg as a function of ratio 
/td is found as shown in Fig. |31 (for details of the 
system parameters see 15j). Similar result has been ob- 
tained using definition (|3b|l of the Ehrenfest time. The 
error bars represent the standard deviations which are 
calculated by taking into account the correlation between 
the first energy levels Ei for systems with slightly differ- 
ent parameters (see0 for details). 
Assuming that 



we found that a = 0.7 ±0.2 and /3 = 0.95 ±0.02. Similar 
results were obtained by using (|3 bp : a — 0.9 ± 0.3 and 
(3 — 1.10±0.07. According to the two competing theories 
the universal values of a = 0.23 and /3 = 1 (Ref. illi) and 
a = 2 and /? = 1 (Ref. Hol) were predicted, while from 
the numerics for Andreev map0 only a = 0.59 ± 0.08 
is universal value. From our numerics the value a is in 
between the two theoretical predictions and within the 
numerical errors it agrees with the result for the Andreev 
map. We found that the values of Ag obtained from (|2b|) 
are slightly greater than those predicted from ljlb|l . How- 
ever, no obvious functional form can be deduced from our 
data for the dependence of Ag on t-e/t^). 

For further support of the effective RMT description 
of the universal gap fluctuation we compare the numer- 
ically obtained ensemble averaged DOS {q{E)) with the 



effective RMT prediction. The ensemble averaged DOS 
{qcs{E)) in the effective RMT description can be calcu- 
lated from (£'rmt(2^)) given in the RMT limit using the 
renormalized parameters Eg and when 'scaling back' 
the variable x ~ {E — Eg)/ Ag into the energy variable E. 

In the RMT limit (g^MTix)) = ~xAi^{x) + [Ai'(x)]^ + 
^ Ai(a::) [l — Ai(y) dy\ is again a universal function 
of X (see note 20 in j6j ) . It can be seen from Fig. 0] that 
the agreement between our numerically obtained ensem- 
ble averaged DOS {g{E)) and {gcs{E)) is excellent at the 
edge of the spectrum (where the theory is valid) for ra- 
tio /td « 0.063 (the same holds for other ratios not 
shown). For larger energies the DOS is around 2/Sn as 
it is expected. For convenience the mean field DOS0 
in the RMT limit is also plotted in the figure. In case 
of two dimensional Andreev billiards the ensemble aver- 
aged DOS {g{E)) shown in Fig. 01 is the first numerical 
evidence for the prediction {gcs{E)) obtained from the 
theoretical result {q-rmt (x)) (see Fig. 2 of Ref. Q). 
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FIG. 4: Ensemble averaged DOS {e{E)) (muhiplied by Sn) 
as a function of E/Et for SA billiard (solid line), the effective 
RMT prediction {qcs{E)) (expressing x in terms of E by us- 
ing Eq. J^J) (short dashes), the Bohr-Sommerfeld calculation 
(dotted fine), and the mean field DOS (long dashes). 

From Fig. ^ one can also see that the agreement be- 
tween the Bohr-Sommerfeld approximation (BS)0 and 
our numerics is quite good at the bottom of the spec- 
trum. This suggests a semi-classical explanation. In a 
work by Schomerus and Beenakker^ close correspon- 
dence has been found between the morphology of the 
phase space and the density of low energy excitations. 
Ref. 4] finds that the semi-classical prediction of the DOS 
for systems with fully chaotic phase space has no definite 
gap but it becomes exponentially small below the energy 
~ 0.5i?T- Furthermore in case of systems with mixed 
phase space and strong coupling of the regular islands to 
the superconductor, the above defined 'gap' is substan- 
tially reduced, namely by a factor Tn/t* , where t* is the 
mean dwell time of trajectories in the chaotic part of the 
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phase space. 

The concept of the gap reduction can also be apphed 
to our SA biUiard, since this is a system with mixed 
phase space (for geometries studied in this work approxi- 
mately 10% of the phase space was regular). To calculate 
the mean dwell time t* in the chaotic part of the phase 
space, only those trajectories are taken into account, for 
which the size of the bunch of trajectories started in its 
Af neighborhood at the N-S interface is increased to the 
characteristic length scale of the billiard before returning 
to the superconducting lead. For the parameters of the 
SA billiard used in our calculations this means that one 
has to exclude trajectories which either do not bounce 
on the circular part or bounce on it only once. From our 
numerics we found (for the ensemble corresponding to 
Fig. 21 that the effective dwell time of these trajectories 
is t* « 1.66td hence the semi-classically obtained gap is 
reduced to the value «0.3i?T- This result is just about 
the energy value where the numerically found DOS be- 
comes exponentially small and thus it confirms the semi- 
classical picture developed in Rcf. 0. Note also that 
below the value « O.SEt (i.e., for x < -2) the RMT 
DOS ( £'rmt(3^)) is also exponentially small. 

While the BS approximation is quite successful in pre- 
dicting the density of low energy excitations, it is clear 
from Fig. ^ that the edge of the spectrum can be better 
predicted using our effective RMT description. This may 
also reveal the limits of the BS approximation. (We ob- 
served similar deviations for other ensembles not shown) . 

In summary we have numerically shown that for finite 
but small enough Ehrenfest time the distribution of the 
rescaled first energy level Ei of an ensemble of chaotic 
Andreev billiards can be treated by an effective RMT de- 
scription. In this model the scaling parameters Eg and 
Ag extracted from the data of the ensemble rescale the 
distribution of Ei such that it agrees with P{x) given in 
the RMT limit. Our numerical results also show that to 
lowest order in t^/tjj the mean field gap Eg decreases 
linearly with the Ehrenfest time but the slope is between 
the two competing theories and close to that of the recent 
numerical calculations for the Andreev map. Calculation 
of the ensemble averaged DOS gives further confirmation 
of our effective RMT description. Our numerics suggest 
that the exponential tail of the DOS can be well inter- 
preted semi-classically. 
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